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1.  Introduction 


Anisotropy  is  prevalent  in  numerous  biological  tissues,  which  are  comprised  of  fibers,  or 
bundles  of  cells  aligned  in  a  uniform  direction.  For  example,  within  the  white  matter 
regions  of  the  human  brain,  axons  tend  to  form  complex  fiber  tracts  that  promote 
anatomical  connectivity  and  communication.  The  bundles  of  axonal  fibers  can  be 
visualized  using  non-invasive  tools  such  as  diffusion  tensor  magnetic  resonance  medical 
imaging  (DTI).  The  basic  principle  behind  DTI  is  that  water  diffuses  more  rapidly  in  the 
direction  aligned  with  the  internal  structure,  and  more  slowly  as  it  moves  perpendicular 
to  the  preferred  direction.  From  the  diffusion  tensor  computed  from  the  imaging, 
diffusion  anisotropy  measures  such  as  the  fractional  anisotropy  (FA)  can  further  be 
computed.  In  addition,  the  principal  direction  of  the  diffusion  tensor  can  be  used  to 
estimate  the  white  matter  connectivity  of  the  brain.  Visualization  of  the  white  matter 
connectivity  is  called  tractography.  An  example  of  DTI  tractography  is  shown  in  figure 
1.  Note  that  DTI  is  one  particular  method  or  application  of  the  broader 
Diffusion- Weighted  Imaging  (DWI). 

The  anisotropic  properties  of  biological  tissue  are  important  to  capture  within  a 
computational  framework  because  the  axonal  fiber  tracts  have  been  reported  to  be 
approximately  three  times  stiffer  than  the  surrounding  matrix  (2).  In  fact,  numerous 
efforts  have  focused  on  incorporating  details  of  fiber  tractography  into  a  computational 
framework  where  many  use  transversely  isotropic  viscoelastic  constitutive  models  to 
represent  the  bulk  and  fiber  mixture.  For  example,  Margulies  et  al.  (2)  conduct 
experiments  and  optimization  techniques  to  develop  a  transversely  isotropic  viscoelastic 
constitutive  model  for  a  porcine  brainstem  tissue  sample.  Wright  et  al.  (3)  present  an 
analytical  and  computational  model  that  also  treats  the  white  matter  as  an  anisotropic, 
hyperelastic  material  using  DTI  to  incorporate  the  structural  orientation  of  the  neural 
axons  into  the  model.  Wright  et  al.  show  that  the  degree  of  injury  that  is  predicted  in  a 
computational  model  of  DAI  is  highly  dependent  on  the  incorporation  of  the  axonal 
orientation  information  and  the  inclusion  of  anisotropy  into  the  constitutive  model  for 
white  matter. 

The  work  presented  here  extends  previous  efforts  by  using  the  transversely  isotropic 
hyperelastic  model  formulation,  and  describes  a  new  design  and  implementation  for 
informing  the  finite  element  model  using  DTI  in  a  three  dimensional,  parallel  computing 
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framework.  Besides  biological  materials,  the  algorithms  described  in  this  paper  can  be 
applied  to  a  broad  range  of  engineering  materials,  as  well.  For  example,  Kevlar  armor 
protection  systems  contain  fiber  directions  that  rotate  with  the  curvature  of  the  surface. 
The  finite  element  implementation  included  here  sets  the  appropriate  framework  for 
incorporating  fiber  directions  into  finite  element  models  used  in  a  diverse  range  of 
applications. 

In  this  report,  we  begin  by  describing  the  continuum  theory  behind  the  implementation 
of  the  transversely  isotropic  hyperelastic  model  in  section  2  followed  by  results  of  a 
validation  study  of  the  implementation  within  the  finite  element  framework  in  section  3. 
Then,  section  4  presents  an  algorithm  for  linking  DTI  and  the  transverse  isotropic  model. 
Section  5  discusses  an  extension  of  the  research  to  modeling  the  white  matter  tissue  of 
the  human  brain  in  a  three-dimensional  finite  element  code.  Then  section  6  offers  some 
concluding  remarks  and  opportunities  for  future  work. 


Figure  1.  Using  Diffusion  Tensor  Imaging  (DTI),  axonal  fiber  tractography  can  be  overlaid  on  a 
finite  element  mesh  and  used  to  construct  a  transversely  isotropic  model.  Tractography  consists 
of  complex  fibers  shown  from  (b)  dorsal,  (c)  right  lateral  side  and  (d)  posterior  views. 
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2.  Continuum  Theory  of  Transversely  Isotropic  Hyperelasticity 


2.1  Kinematics  of  Finite  Deformations 


We  begin  by  defining  the  deformation  gradient  defined  as: 


(1) 


where  x  is  the  deformed  configuration  of  a  material  particle  and  X  is  the  reference 
position.  The  ratio  of  deformed  to  undeformed  volumes  is  given  by  the  determinant  of 

F, 

J  =  det(F)  (2) 

Some  materials  behave  differently  in  bulk  and  shear,  so  it  is  most  beneficial  to  perform  a 
multiplicative  decomposition  of  F  into  volume-changing  (dilational)  and 
volume-preserving  (distortional)  parts  (4).  To  do  this,  a  deviatoric  deformation  gradient 
is  defined  in  which  the  volume  change  is  eliminated: 


F  =  J~^3F 


(3) 


It  then  follows  that: 

C  =  FTF  (4) 

B  =  FFT  (5) 

where  C  and  B  are  referred  to  as  the  modified  right  Cauchy-Green  tensor,  and  the 
modified  left  Cauchy-Green  tensor,  respectively  The  modified  principle  invariants  of 
the  Cauchy-Green  deformation  tensor  are: 

7,  =  tr  C  =  tr  B 

I2  =  i  [(tr  C)2  -  tr  (C2)]  =  i  [(tr  B)2  -  tr  (b2)[ 

J3  =  det  C  =  det  B 


(6) 

(7) 

(8) 
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Two  additional  invariants  are  introduced  in  order  to  capture  the  mechanics  of  transverse 
isotropic  materials: 


I  a  =  a0  ■  Ca0  (9) 

I5  =  a0-C2a0  (10) 

where  J4  and  I5  arise  directly  from  the  anisotropy  and  describe  the  properties  of  the 
fiber  family  and  its  interaction  with  the  other  material  constituents.  The  vector  a0  is 
introduced  to  describe  the  undeformed  fiber  direction  where  we  have  assumed  that  the 
only  anisotropic  property  of  the  solid  comes  from  the  presence  of  the  fibers.  For  a 
material  that  is  reinforced  by  only  one  family  of  fibers,  the  stress  at  a  material  point 
depends  not  only  on  the  deformation  gradient  F,  but  also  on  that  single  preferred 
direction,  which  is  called  the  fiber  direction. 

2.2  Decoupled  Representation  of  Strain  Energy 

Following  the  approach  taken  by  others  (2,  4,  5)  in  the  representation  of  transversely 
isotropic  hyperelasticity,  a  unique  decoupled  representation  of  the  strain  energy  function 
per  unit  volume  in  terms  of  the  five  invariants  is  given  by: 


^(Ili  1 2 j  ^3,  I4,  1 5)  —  'Tiso(/1,  J2,  1 3)  +  Tanjso (/4,  1 5) 


(11) 


where  T'iso(J 4,  J2,  / 3 )  describes  the  response  of  the  isotropic  matrix,  and  T;iniS0(I4.  h) 
describes  the  directional  contribution  of  the  fiber  reinforcement.  Note  that  for  an 
incompressible  material,  J3  =  J2  =  1. 


2.3  Uncoupled  Stress  Tensor 


For  a  hyperelastic  material,  the  2nd  Piola-Kirchhoff  stress  is  derived  from  the  strain 
energy  as: 


S 


2 


<9T 

dC 


2J-2/3DEV 


<9T  (C) 
dC 


-pj  C1 


(12) 

(13) 


where  the  operator  DEV  [•]  =  [•]  —  |  ([■]  :  C)  C  1  extracts  the  deviatoric  part  of  a  tensor  in 
the  reference  configuration  and  p  =  —  is  the  hydrostatic  pressure.  Using  the  chain 
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rule,  equation  13  can  be  further  written  as: 


S  =  2J-2/3DEV 


-pj  C”1 


(14) 


where  the  terms  jh  are  given  as: 

(15) 

(16) 

(17) 

(18) 

where  I  is  a  2nd  rank  identity  tensor. 


dh 

dC 

dl2 

W 

dl4 

w 

dh 

oc 


^4  =  7,1  -  C 


=  =  a0  (g)  a0 


-z=  =  a0  <g>  C  ■  a0  +  a0  •  C  (8)  a0 


A  push-forward  operation  on  the  second  Piola-Kirchoff  stress  tensor  S  to  the  current 
configuration,  and  a  scaling  with  the  inverse  of  the  volume  ratio,  transforms  equation  16 
to  the  Cauchy  stress  tensor  of  a  transversely  isotropic,  compressible  hyperelastic 
material: 


a  =  —dev 

U 


=<94'  —  t 

F^F 

dC 


—  pi 


(19) 


a  =  —dev 
u 


^4q  -/34c)  8  —  'k2B"J  -|-  /44>4a  (^)  a  744' 5  ^a  <£<)  Ba  -)-  aB  <£<)  a) 


pi  (20) 


where 


dev  [•] 

% 

a 


[•]  -  7.  (H  :  I)  I 

dh 

Fa0 


(21) 

(22) 

(23) 


in  which  a  describes  the  configuration  of  a0  during  the  deformation. 
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2.4  Uncoupled  Stress  Tensor  for  a  Particular 

Within  the  framework  of  incompressible,  transversely  isotropic  hyperelasticity,  the  strain 
energy  function  can  be  written  in  the  uncoupled  form  as: 

*(C)  =  *voi(J)+^(C)  (24) 

where  'EVoi(  J)  represents  the  volumetric  strain  energy  contribution  and  4qsoc  (C) 
represents  the  isochoric  strain  energy  function.  Following  Weiss  (5),  the  isochoric  part 
of  the  strain  energy  can  be  further  subdivided  into  separate  functions  associated  with 
the  modified  invariants: 


tf(C)  =  ^voi(d)  +  F 1  (J1;  I2)  +  F-2  (J4)  (25) 

with: 

F \  =  Ci  (h  -  3)  +  U2  (J2  -  3)  (26) 

F2  =  C3  [exp  (J4  -  1)  -  I4]  (27) 


where  C\,  C2,  and  6'3  are  constants  obtained  from  parameter  fits  to  experimental  data. 
Equation  26  corresponds  to  a  modified  Mooney-Rivlin  model,  and  equation  27  describes 
an  exponential  behavior  in  the  fiber  direction  that  is  a  characteristic  of  most  soft  tissues. 


Then,  similar  to  equation  13,  the  second  Piola-Kirchoff  stress  is  given  by: 


S 


2J-2/3DEV 


^iSQC  (C) 
dC 


d^Voi  (J) 

dJ 


j  cr1 


(28) 


where  the  partial  derivative  of  the  isochoric  part  is  given  by: 


d^isoc  (c) 


dC 


d^isoc  (C)  dtf  isoc  (C)7 
H - ^ - J-  i 


dh 


dl2 


j  dd'isoc  (C)  —  d^isoc(C) 

I - FF~  ~C  + - Fi~ — “a°  ®  a°  (29> 

Ol  2  Ul 4 
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so  then  substituting  equation  29  into  the  DEV  operator  in  equation  28, 


DEV 

<9^isoc  (C) 

d^isoc{C)  d'hjsoc  (C)  - 

—  —  1  1 

dc 

dh  dl2 

<9 'Eis 


<9  I 


(clu+^iS0C  (C) 


dh 


a0  (8)  a0 


1 

3 


<9 'his 


dh 


(ch1+2d*iaoc 


dl 


Mi2  + 


<9 'hjs 


dl. 


(c)/ 

- -*4 


c 


-1 


(30) 


Then,  the  second  Piola-Kirchoff  stress  is  given  by: 


S  =  2  J~2/3 


'd*isoc  (C)  d'h.soc  (C)  \  d^SOC  (C)  d'hisoc  (C) 

dh  dl  2  j  dh  dl4 


1  ( <9Vhisoc  (C)  +  o^'hisoc  (^)  J  ^  _|_  ^isoc  (C)  |  ^_,_1 


dh 


dh 


dh 


«9^vol  (J) 
dJ 


a0  8)  a0 


J  cr1 


(31) 


where. 


d'hvo!  ( J) 

dJ 

d%^c{C) 

dh 

<9^isoc  (C) 
dl2 

d%hc  (C) 

dl4 


V 

Cl 

c2 

C3  [exp  (74  -  l)  -  l] 


(32) 

(33) 

(34) 

(35) 


Then,  substituting  into  equation  31: 

S  =  -pJC~1+ 

2  J"2/3  [(Ci  +  Cali)  I  -  C2C  +  C3  [exp  (I4  -  l)  -  l]  a0  8  a0  (36) 

-1/3  (<7,7,  +  2C2/2  +  C3  [exp  (74  -  1)  -  1]  74)  C-1] 


a  push-forward  operation  yields  the  Cauchy  stress: 


a  =  —pl+ 


2 

J 


{Ci  +  C2h)  B 


C2B~  +  C3  [exp  [h  -  l)  -  l]  a  8  a 


-1/3  {Cih  +  2 C2/2  +  C3  [exp  {IA  -  1)  -  1]  h)  I 


(37) 
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3.  Validation  of  Numerical  Method 


To  insure  that  the  stress  update  was  properly  implemented  into  the  computational 
framework,  numerical  tests  were  performed  and  compared  to  an  analytical  solution. 
Single  element  computations  of  uniaxial,  strip  biaxial,  and  equibiaxial  tests  were 
performed  to  validate  the  update  of  the  stress  for  the  implementation  against  theoretical 
solutions.  If  the  stresses  are  applied  in  the  x-y  plane  along  the  principal  axes  so  the  azz 
component  is  zero  and  the  incompressibility  constraint  det  F  =  1  is  applied,  the 
deformation  gradient  can  be  expressed  as: 


F 


\x  0  0 

0  Xy  0 


(38) 


Then  using  equation  37,  the  stresses  are: 


(Cl  +  C2h)  (Bxx  -  B„)  -  C2  (bL  -  Bt)  + 


C3  [exp  (J4  -  l)  -  l]  1 4  (a2x  -  a2z) 


=  2  [(C!  +  C2h)  (B„  -  B„)  -  C2  (b2„  -  B]z)  + 
C3  [exp  (74  -  l)  -  l]  J4  (a2  -  a2z)\ 


(39) 


(40) 


In  the  case  of  uniaxial  extension,  ayy  =  azz  =  0  and  Xy  =  Xz.  For  a  strip  biaxial  test, 
azz  =  0  and  A2  =  1.  For  an  equibiaxial  test,  azz  =  0  and  Xx  =  Xy.  In  addition,  the  same 
material  constants  were  used  for  all  three  cases  including,  C\  =  10.0,  C2  =  10.0, 

C's  =  100.0  and  n  =  10.0e9.  As  seen  in  figure  2,  good  agreement  between  the  theoretical 
and  finite  element  solutions  for  all  three  cases  of  loading  is  acheived.  The  convergence 
of  the  finite  element  solution  to  the  theoretical  solution  demonstrates  that  the  stress 
update  equations  were  implemented  correctly. 
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Figure  2.  Comparison  between  the  finite  element  results  and  analytic  solutions  for  three  different 
loading  conditions,  (a)  uniaxial,  (b)  strip  biaxial  and  (c)  equibiaxial. 


4.  An  Algorithm  for  Linking  DTI  and  Continuum  Mechanics 


The  transversely  isotropic  hyperelastic  model  previously  implemented  depends  on  the 
fiber  reinforcement  in  the  undeformed  configuration,  a0.  The  algorithm  developed  here 
describes  how  to  assign  each  finite  element  an  orientation  based  on  the  fiber 
tractography  obtained  from  DTI,  as  introduced  in  section  1.  This  is  complex  because  of 
the  non-linear  tractography  and  unstructured  finite  element  shapes  within  the 
three-dimensional  framework.  Algorithms  for  both  hexahedral  and  tetrahedral 
elements  are  described. 

First,  the  algorithm  divides  each  fiber  line,  /*  (where  i  goes  between  0  and  the  number  of 
tractography  fibers),  into  individual  fiber  segments,  s j  (where  j  goes  between  1  and  the 
number  of  segments  within  a  given  fiber  line).  Each  fiber  segment  consists  of  two  fiber 
segment  points  with  coordinates  pi  and  p2:  s  =  p2  —  p\.  The  first  discretized  point 
within  the  fiber  segment  is  assigned  a  direction  based  on  the  vector  connecting  it  to  the 
second  point.  If  the  fiber  point  exists  at  the  end  of  a  fiber  line,  that  terminating  point 
receives  the  same  orientation  as  the  previous  point.  A  schematic  of  the  discretization 
process  is  shown  in  figure  3. 

Next,  the  algorithm  determines  if  a  fiber  segment  within  a  given  tractography  fiber 
overlaps  spatially  with  a  finite  element.  In  order  to  determine  this,  various  cases  need  to 
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Fiber  point  receives 


Figure  3.  Fiber  discretization  and  notation  definitions  used  throughout  the  report. 


be  analyzed.  We  begin  our  search  algorithm  at  the  fiber  segment  end  points.  For  each 
element  in  the  volumetric  mesh,  every  point  within  each  fiber  is  analyzed  for  the  cases 
schematically  shown  in  figure  4.  Case  1  exists  if  the  first  point  in  a  fiber  segment  is 
contained  in  a  given  element.  Case  2  exists  if  the  second  point  of  a  fiber  segment  is 
contained  in  a  given  element.  Case  3  exists  if  both  endpoints  of  the  fiber  segment  are 
contained  within  a  given  element.  Case  4  exists  if  a  fiber  segment  intersects  the  faces  of 
the  element,  but  its  fiber  points  are  not  contained  within  a  given  element.  Case  5  exists 
if  the  element  contains  none  of  the  fiber  points,  and  the  segment  does  not  intersect  the 
element  faces.  Note,  that  there  is  an  unlikely  subset  of  case  3  when  a  given  fiber 
segment  lies  parallel  with  a  single  element  face,  so  this  orientation  would  be  shared  by 
any  adjacent  element  also. 

For  all  of  the  cases  shown  in  figure  4,  determining  if  a  fiber  point  is  contained  within  a 
finite  element  is  complicated  by  the  unstructured  mesh  shapes.  For  example,  in 
elements  with  adjacent  faces  far  from  orthogonal  with  each  other,  a  fiber  point  may  lie 
outside  of  the  faces  of  the  element  even  if  the  fiber  point  spatial  coordinates  are  in  the 
bounds  of  the  maximum  and  minimum  element  nodal  coordinates.  Therefore,  an 
algorithm  to  determine  if  the  fiber  point  is  inside  of  the  element  is  implemented.  The 
algorithm  takes  into  consideration  the  angles  of  the  element  faces  explicitly  by 
computing  normals  to  the  element  faces  and  the  coordinates  of  the  fiber  point,  as 
schematically  shown  in  figure  5. 

A  slightly  more  complicated  algorithm  is  used  for  case  4  where  the  line  segment 
intersects  the  element  faces,  but  no  points  in  the  fiber  segment  reside  within  the  element. 
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Figure  4.  Schematic  of  the  various  cases  that  are  possible  when  determining  if  a  fiber  segment 
(within  a  given  tractography  fiber)  overlaps  with  a  finite  element  in  space. 


Tetrahedron  Elements 


2 


A  x  B  *D  AND  A  x  B  «P 
F  x  E  •  B  AND  F  x  E  •  P 
G  xD  •EANDGxD'P 
BxD  •A  AND  BxD.p 


If  each  pair  has  the  same  sign, 
then  the  point  is  inside  the 
tetrahedron 


Hexahedron  Elements 


S’x'B  *DANDXxB  •'F 

FxTi  *Ti andT x_h_*p 

^_X  K  *H  AND  I_X  KjP 
LxD-KANDrxD  *T 
FxM'lANDlxM'P 
TT  x"B  »T  ANDTT  x^.T 


If  each  pair  has  the  same  sign,  then 
the  point  is  inside  the  hexahedron 


Figure  5.  Method  used  to  test  if  a  specified  point  is  within  the  element.  Algorithms  for  both 
tetrahedron  and  hexahedron  elements  are  designed  and  implemented. 
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As  schematically  shown  in  figure  6,  first  the  algorithm  calculates  the  normal  of  the  first 
face,  n  by  taking  the  cross  product  of  two  vectors  which  share  the  same  vertex.  Pi. 

n  =  (P-2  —  Pi)  x  (P:i  -  Pi)  (41) 

Next,  the  algorithm  takes  the  dot  product  of  the  distance  between  an  endpoint,  Li  or  L2, 
and  the  shared  vertex  Pi.  This  results  in  two  dot  products,  the  distances  from  Li  and  L> 
from  P] : 

xi  =  (Li  -  Pi)  •  (n)  (42) 

%2  =  (L2  -  Pi)  •  (n)  (43) 

If  xi  ■  x2  >  0,  the  line  segment  never  intersects  the  plane  of  the  face.  If  x\  =  x2,  the  line 
segment  and  the  plane  are  parallel,  and  the  line  segment  can  be  disregarded  since  it  does 
not  intersect  two  different  faces  of  the  element.  If  the  two  preceding  conditions  do  not 
exist,  the  code  will  find  the  location  of  the  point,  P,  on  the  given  line  segment  which 
intersects  the  plane. 

P  =  L1  +  (L2-L1)(^~)  (44) 

x2  -  Xi 

The  point  P  is  the  intersection  point  that  exists  on  the  plane  of  the  element  face;  therefore, 
the  algorithm  must  test  if  P  also  lies  within  the  surface  area  of  an  element  face  (see 
figure  6).  A  point  is  inside  of  a  polygon  if  it  is  always  on  the  same  side  of  all  the 
segments  that  create  a  continuous  loop  around  the  edges.  For  a  given  fiber  segment  that 
solely  intersects  the  element  faces,  there  should  be  exactly  two  intersection  points  on  a 
single  element. 

If  there  exist  elements  for  a  given  finite  element  mesh  where  no  fiber  tracks  overlap,  the 
elements  are  assigned  an  orientation  vector  of  zero,  a0  =  0.  Therefore,  the  transversely 
isotropic  material  condenses  to  an  isotropic  material. 

4.1  A  Scheme  for  Averaging  Multiple  Fibers 

The  previous  discussion  described  the  methods  used  to  find  the  existence  and  direction 
of  a  single  fiber  in  a  hexahedral  or  tetrahedral  element.  The  fiber  direction  assigned  to 
an  element  with  only  one  fiber  segment  is  given  by: 

a°  =  TTTT  (45) 

s 
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Figure  6.  Method  used  to  test  if  a  line  intersects  the  faces  of  an  element  for  the  case  in  which  the 
line  segment  intersects  the  element  faces,  but  no  points  in  the  segment  reside  within  the  element 
(See  Case  4  in  figure  4). 
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where  s  represents  the  orientation  vector  of  the  fiber  segment.  Depending  on  the  fiber 
or  mesh  densities,  there  could  be  multiple  fibers,  multiple  fiber  segments,  or  both 
existing  in  a  single  element,  as  schematically  depicted  in  figure  7.  The  current  effort 
develops  an  average  scheme  for  each  possibility  since  the  current  implementation  only 
permits  one  orientation  vector  per  element.  In  future  work,  we  hope  to  extend  the 
model  to  have  the  capability  to  include  multiple  fiber  orientations  within  a  single  finite 
element,  which  will  require  additional  modifications  to  the  theoretical  basis.  The 
method  of  Diffusion  Spectrum  Imaging  (DSI)  is  an  example  of  where  this  functionality 
would  be  useful  to  have,  since  DSI  is  capable  of  describing  fiber  crossings. 
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Element  Intersection 


71 
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Figure  7.  Depending  on  the  fiber  or  mesh  densities,  there  could  be  multiple  fibers,  or  multiple 
fiber  segments,  or  both  existing  in  a  single  element  which  requires  an  averaging  scheme. 


In  the  case  of  a  single  fiber  with  multiple  fiber  segments  within  the  same  element,  an 
approximation  of  an  "average"  vector  direction  is  taken  because  the  material  model  uses 
only  one  fiber  direction  per  element.  By  simply  adding  the  vectors,  you  obtain  a 
resultant  "average"  direction,  which  provides  an  estimate  of  an  overall  fiber  orientation 
as  schematically  depicted  in  figure  8. 

For  two  or  more  fibers  in  a  single  element,  the  assigned  orientation,  a0,  is  calculated  by 
the  following: 
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max 

max  segments 
fibers  in  fiber  i 


sbiw 


a, 


(46) 


p2 


P3 


Figure  8.  The  sum  of  the  vectors  will  give  an  average  orientation. 


Uf)  — 


(47) 


where  s j  are  the  orientation  vectors  of  the  fiber  segments  within  the  element,  and  ae 
represents  the  average  of  s 3  -  two  averages  are  taken  with  each  element.  First,  the 
vectors  from  the  segments,  s3,  within  the  same  fiber  line  and  element  are  averaged. 
Second,  the  average  of  all  fibers  within  that  element  are  calculated.  Because  the  fiber 
average  scheme  depends  on  the  element  size,  future  efforts  will  explore  the  effects  of 
finite  element  mesh  density.  A  course  finite  mesh  decreases  computational  cost,  but 
may  also  smear  out  the  complex  fiber  orientations.  However,  there  are  regions  of  the 
fiber  tractography  that  change  slowly  over  spacial  distance  so  a  course  mesh  may  be 
more  appropriate.  The  preceding  procedures  for  determining  an  element  orientation  for 
the  transverse  isotropic  model  can  be  summarized  by  the  global  algorithm  shown  in 
figure  9. 


for  each  element 
for  each  fiber 

for  each  segment  in  fiber 

Determine  case  for  fiber  points 
Store  segment  orientation  vector  s 
end 

Sum  segment  vectors  for  given  element  and  fiber 

end 

Compute  ae  (Sum  orientation  vectors  for  all  fibers  in  element) 
Compute  a0 


Figure  9.  The  global  algorithm  used  to  obtain  final  vector  orientations  for  each  element. 
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Although  this  method  will  automatically  take  into  account  the  vector  magnitudes,  one 
limitation  of  this  approach  involves  the  loss  of  orientation  information  due  to  the 
averaging  scheme  for  adjacent  fiber  segments  that  are  close  to  perpendicular.  For 
example,  as  shown  in  figure  10,  if  one  point  of  the  fiber  segment,  plf  is  entirely  outside  of 
the  element,  while  the  other  point  of  the  segment,  p2,  is  inside  the  element  but  relatively 
close  to  the  element's  surface,  that  fiber  segment  direction  will  receive  more  weight  than 
it  should.  It  receives  the  same  amount  of  weight  as  it  would  if  it  were  fully  contained 
within  the  element.  Although  the  accuracy  of  the  fiber  representation  is  limited  by  the 
mesh  resolution,  it  is  possible  in  the  future  to  create  a  more  accurate  averaging  of  fiber 
directions  in  each  element. 


Pi 


Figure  10.  Averaging  fibers  which  are  not  fully  contained  in  the  element  can  create  a  misrepre¬ 
sented  average  direction  for  the  individual  element. 


To  visualize  the  assigned  element  fiber  directions,  the  unit  vector  was  scaled  by  a 
characteristic  length  of  the  element  (diagonal  or  side  length)  and  added  to  the  center 
coordinates  of  each  element.  The  VTK  output  file  describes  the  lines  connecting  the 
center  positions  and  endpoints  of  the  scaled  unit  vectors,  which  can  be  opened  up  in 
Paraview  (Kitware)  and  then  compared  against  the  original  dataset.  The  visualization 
lines  are  not  meant  to  completely  coincide  with  the  input  fiber  lines  but  rather  represent 
the  directions  assigned  to  each  element.  Since  the  visualization  lines  are  positioned  at 
the  center  of  each  element,  the  lines  will  be  in  close  proximity  to  the  fiber  lines  they 
represent,  but  will  not  be  exactly  overlapping.  A  higher  resolution  mesh  will  produce  a 
more  accurate  representation  of  fiber  direction  assignment  to  the  elements. 
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5.  Application  to  Modeling  White  Matter  Tissue  of  the  Human  Brain 


The  primary  objective  of  implementing  the  transverse  isotropic  model  was  to 
incorporate  fiber  directions  into  the  finite  element  models  for  application  to  anisotropic 
biological  tissue.  Multiple  scenarios  were  created  to  test  the  algorithms;  however, 
instead  of  presenting  the  exhaustive  list  here,  only  the  brain  model  is  shown.  Figure  11a 
shows  the  input  fiber  tractography  superimposed  with  a  finite  element  mesh  of  the 
human  brain.  Figure  lib  shows  the  fiber  orientations  assigned  to  the  elements,  obtained 
from  the  fiber  tractography.  While  this  example  shows  qualitative  agreement,  it  is 
important  to  note  that  the  finite  element  mesh  and  the  fiber  tractography  are  not  from 
the  same  individual  at  this  time.  There  is  currently  a  major  focus  to  develop  models 
from  a  single  individual. 
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6.  Concluding  Remarks 


The  work  presented  here  provides  a  novel  way  to  use  state-of-the-art  medical  imaging  to 
inform  a  3-D  transversely  isotropic  model  of  biological  tissue.  After  testing  the 
numerical  model's  implementation,  an  algorithm  was  created  in  order  to  read  fiber  data 
and  output  the  undeformed  fiber  directions  into  the  material  model.  The  initial  results 
are  promising  since  they  qualitatively  show  adequate  correlation  between  the  original 
DTI  data  and  the  averaged  directions  that  were  assigned  to  the  elements.  The 
DWI-informed  modeling  can  be  extended  to  other  biological  tissues  and  different 
material  models.  Aside  from  using  DTI  for  imaging  the  white  matter  of  the  brain,  it  can 
also  measure  skeletal  muscle  fiber  directions  in  fixated  skeletal  muscles.  Additionally, 
the  model  could  be  useful  for  finite  element  modeling  of  non-biological  fibrous 
materials,  such  as  kevlar  and  various  synthetic  polymers. 

As  mentioned  earlier,  future  work  will  include  finite  element  models  that  correspond  to 
the  MRI  data  and  DTI  data  derived  from  the  same  individual.  The  algorithm  could  be 
improved  by  constantly  updating  the  fiber  directions,  which  may  change  during 
deformation,  although  this  would  most  likely  increase  the  overall  time  required  to  run 
the  simulation  and  may  not  be  necessary  for  simulations  of  short  events,  such  as  blast  or 
impact  loading.  Furthermore,  DSI  images  may  be  incorporated  into  the  computational 
framework  so  that  fiber  crossings  may  also  be  included. 
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UNIVERSITY  OF  NEBRASKA 
114G  OTHMER  HALL 
PO.  BOX  880642 
LINCOLN,  NE  68588-0642 
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1  LAKIESHA  N.  WILLIAMS 

MISSISSIPPI  STATE  UNIVERSITY 
BOX  9632 

MISSISSIPPI  STATE,  MS  39762 

1  BARCLAY  MORRISON 
COLUMBIA  UNIVERSITY 
351  ENGINEERING  TERRACE 
1210  AMSTERDAM  AVENUE,  MAIL 
CODE:  8904 
NEW  YORK,  NY  10027 

1  SHANE  SCHUMACHER 

SANDIA  NATIONAL  LABORATORIES 

NANOSCALE  AND  REACTIVE 

PROCESSES 

PO.  BOX  5800,  MS  0836 

ALBUQUERQUE,  NEW  MEXICO 

87185-0836 

1  CHARLES  E.  NEEDHAM 
APPLIED  RESEARCH  ASSOCIATES 
INC. 

SOUTHWEST  DIVISION 
4300  SAN  MATEO  BLVD.,  NE 
SUITE  A-220 

ALBUQUERQUE,  NM  87110 

2  ARMY  RSRCH  OFC 
RDRL  ROE  M 

D  STEPP 
RDRL  ROE  V 
L  RUSSELL 
PO  BOX  12211 

RESEARCH  TRIANGLE  PARK,  NC 
27709-2211 

1  ADAM  FOURNIER 

U.S.  ARMY  ABERDEEN  TEST 
CENTER 

ATTN:  TEDT-AT-SLB 
400  COLLERAN  ROAD 
ABERDEEN  PROVING  GROUND, 
MARYLAND  21005-5059 

1  UNIV  OF  FLORIDA 

MECHL  AND  AEROSPACE  ENGRNG 
G  SUBHASH 
GAINESVILLE,  FL  32611 
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5  THE  JOHNS  HOPKINS  UNIVERSITY 
APPLIED  PHYSICS  LABORATORY 
11100  JOHNS  HOPKINS  ROAD 
LAUREL,  MARYLAND  20723-6099 
ANDREW  MERKLE 
MORGANA  TREXLER 
ANDY  LENNON 
JACK  ROBERTS 
TIM  HARRIGAN 

2  MASSACHUSETTS  INSTITUTE  OF 
TECHNOLOGY 
INSTITUTE  FOR  SOLDIER 
NANOTECHNOLOGIES 
BLDG  NE47, 4TH  FLOOR 
77  MASSACHUSETTS  AVENUE 
CAMBRIDGE,  MA  02139 
RAUL  RADOVITZKY 
SIMONA  SOCRATE 

1  PAUL  E.  RAPP 

DIRECTOR,  TRAUMATIC  INJURY 
RESEARCH  PROGRAM 
DEPARTMENT  OF  MILITARY  AND 
EMERGENCY  MEDICINE 
UNIFORMED  SERVICES  UNIVERSITY 
OF  THE  HEALTH  SCIENCES 
4301  JONES  BRIDGE  ROAD 
BETHESDA,  MARYLAND  20814-4799 

1  JOHN  M.  GETZ 

U.S.  ARMED  FORCES  MEDICAL 
EXAMINER  SYSTEM 
1413  RESEARCH  BLVD. 

ROCKVILLE,  MARYLAND  20850 

1  RUTGERS  THE  STATE  UNIV 
OF  NEW  JERSEY 
DIR  OF  FED  RSRCH  RELATIONS 
B  LAMATTINA 
96  FRELINGHUYSEN  RD 
CORE  BLDG 
PISCATAWAY  NJ  08854 

1  ALAN  HEPPER 

DSTL  BIOMEDICAL  SCIENCES 
RM  1A,  BLDG.  245 
PORTON  DOWN 
SALISBURY,  WILTSHIRE 
SP4  OJQ 
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4  DRDC  VALCARTIER 

2459,  PIE-XI  BLVD.  NORTH 
QUEBEC,  QC  G3J  1X5  CANADA 
KEVIN  WILLIAMS 
AMAL  BOUAMOUL 
LUCIE  MARTINEAU 
DENNIS  NANDLALL 

1  CRAIG  BURRELL 
DRDC  TORONTO 
1133  SHEPPARD  AVENUE  WEST 
PO.  BOX  2000 

TORONTO,  ON  M3M  3B9,  CANADA 

1  ALICE  B.  AIKEN 
CIMVHR 

QUEEN'S  UNIVERSITY 
SCHOOL  OF  REHABILITATION 
THERAPY 

KINGSTON,  ON  K7L3N6,  CANADA 

2  CENTER  FOR  INJURY 
BIOMECHANICS 

WAKE  FOREST  UNIVERSITY 
MEDICAL  CENTER  BLVD. 
WINSTON-SALAM,  NC  27157 
JOEL  STITZEL 
F.  SCOTT  GAYZIK 

1  DANIEL  WISE 

HENRY  JACKSON  FOUNDATION 
U.S.  ARMY  AEROMEDICAL 
RESEARCH  LABORATORY 
6901  ANDREWS  AVENUE 
FORT  RUCKER,  AL  36362-0577 

1  TOM  RADTKE 

HUMAN  PROTECTION  AND 
PERFORMANCE  DIVISION 
DEFENCE  SCIENCE  AND 
TECHNOLOGY  ORGANISATION 
DEPARTMENT  OF  DEFENCE 
BLDG  109,  506  LORIMER  STREET 
FISHERMANS  BEND,  VICTORIA  3207 
AUSTRALIA 
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2  DEPARTMENT  OF  MECHANICAL 
ENGINEERING 

THE  JOHNS  HOPKINS  UNIVERSITY 
LATROBE  122 

3400  NORTH  CHARLES  STREET 
BALTIMORE,  MD  21218 
K.T.  RAMESH 
VICKY  NGUYEN 

1  GARY  KAMIMORI 

DEPARTMENT  OF  BEHAVIORAL 
BIOLOGY 

WALTER  REED  ARMY  INSTITUTE  OF 
RESEARCH 

503  ROBERT  GRANT  AVENUE,  2W97 
SILVER  SPRING,  MD  20910-7500 
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